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This paper offers a reappraisal of Fung's model 
for quasi-linear viscoelasticity. It is shown that a 
number of negative features exhibited in other works, 
commonly attributed to the Fung approach, are 
merely a consequence of the way it has been applied. 
The approach outlined herein is shown to yield 
improved behaviour and offers a straightforward 
scheme for solving a wide range of models. Results 
from the new model are contrasted with those in 
the literature for the case of uniaxial elongation of 
a bar: for an imposed stretch of an incompressible 
bar and for an imposed load. In the latter case, a 
numerical solution to a Volterra integral equation is 
required to obtain the results. This is achieved by a 
high-order discretization scheme. Finally, the stretch 
of a compressible viscoelastic bar is determined for 
two distinct materials: Horgan-Murphy and Gent. 

1. Introduction 

The study of nonlinear viscoelastic deformations of solid 
materials has a very long history, with a consequent 
proliferation of a diverse and extensive array of 
constitutive models. For a comprehensive overview 
of the topic, the reader is directed to the recent paper 
by Wineman [1] and related references therein. The 
subject is still active, with models continuing to be 
developed across the field, from highly mathematical 
approaches where implementation is not a concern 
to very applied studies where ease of application 
is essential. Each approach has its advantages: the 
models that more accurately describe the microphysics 
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tend to prove difficult to employ in practical engineering or biomedical situations, whereas 
simpler approaches often miss crucial details. If the viscous model assumes a fading memory 
effect of the strain history, then this usually gives rise mathematically to Volterra-type integral 
equations, which in the linear case can generally be solved either analytically or numerically 
However, in the nonlinear case (of interest to describe finite deformations), these problems are 
difficult to solve by any methods, even finite-element approaches. It is therefore crucial to develop 
constitutive models that are simple enough to be amenable to straightforward (and rapid) 
solution methods, yet include enough detail to capture the important physics underlying these 
relevant materials. One such 'compromise' approach was offered by Fung [2], in order to study 
the uniaxial elongation of biological soft tissue. His constitutive assumption, often called quasi- 
linear viscoelasticity (QLV) or Fung's model of viscoelasticity, assumes that the viscous relaxation 
rate is independent of the instantaneous local strain. This model, which is a special case of a more 
general Pipkin-Rogers constitutive model [3], predicts that at any time a stress that is equal to 
the instantaneous elastic stress response decreased by an amount depending on the past history, 
assuming that a Boltzmann superposition principle holds. 

Fung's QLV model is perhaps the most widely used today, but has met with some criticism 
mainly because it has been suggested that it does not always yield physically 'reasonable' 
behaviour. Of course, as with any constitutive model for a complex nonlinear material, QLV has 
limitations. However, it can be expected to be appropriate for materials whose relaxation-rate 
coefficients are weakly dependent on deformation, or where the deformation comprises small 
perturbations about a large initial deformation. Various experimental results have appeared in 
the literature that appear to confirm that QLV is, in practice, a reasonable model for a range of 
materials; see for example [4-6]. 

The apparent deficiency of QLV will be discussed below, but the purpose of this paper 
is to reappraise Fung's approach, suggesting how it can be reformulated against existing 
interpretations in the literature. It will be seen that the form offered herein has similarities to other 
viscoelastic formulations, including that by Simo [7], and most importantly exhibits behaviour 
that makes it both (relatively) easy to use and physically useful. 

Analysis will commence with the general tensorial form of (QLV) proposed by Fung [2] 



n(x,o = 



G(f-r) dr, (1.1) 



where G{t) is, in comparison with linear theory, the stress relaxation second-order tensor, TT 
denotes the second Piola-Kirchhoff stress, while TT^ is an instantaneous strain measure. The latter 
may be thought of as an equivalent (instantaneous) elastic stress, hence the superscript 'e'. This 
tensorial integral identity is the natural generalization of the simple one-dimensional relationship 
proposed by Fung [2], which preserves objectivity. 

As mentioned, thanks to the relative simplicity of Fung's approach over more general 
nonlinear viscoelastic models, it has proved extremely popular. This is especially true in the 
biomechanics field where it has been employed to predict the large deformations of soft tissues. 
However, despite the widespread use of QLV, published studies (especially in the incompressible 
limit) appear to interpret (1.1) in different ways, which then lead to quite different results. In the 
paragraphs below, a number of common approaches are discussed, where the interpretations of 
the model may be questioned. The main aim of this paper is then to re-derive QLV, starting from 
basic principles and also to ensure consistency in the limit of infinitesimal deformations, thus 
recovering the Boltzmann theory of linear viscoelasticity. 

The simplest interpretation of Fung's relation is to assume a purely one-dimensional 
homogeneous deformation. This is attractive for its ease of use and has been employed 
extensively in biomedical applications [8-10], for example to model the shearing deformation 
of brain tissue [6]. However, in practical applications, especially for incompressible or near 
incompressible materials purely one-dimensional deformations, such as pure uniaxial extension, 
will not be realizable. Thus, a tensorial form of QLV must be employed. As mentioned. 



this relation must, at the minimum, satisfy objectivity, and hence the form chosen in (1.1). 
Nevertheless, a number of authors have erroneously expressed Fung's relation not for the second 
Piola-Kirchhoff stress, which guarantees objectivity, but in terms of other stress measures. For 
example, see Fung's original discussion on the subject in [2, page 253], in which the QLV relation 
is expressed in terms of the Kirchhoff stress! In the latest version [6.12] of ABAQUS, although the 
viscoelastic constitutive model is objective, the model is rather heuristic and as we shall see later, 
it assumes the same relaxation behaviour in both the shear and bulk parts of an isotropic material. 
Further technical aspects of objectivity are discussed in the paper by Liu [11]. 

Another approach to QLV sometimes employed in the literature (e.g. [12-14]) is, in an 
analogous fashion to nonlinear elasticity, to write the stress in terms of the instantaneous 
derivative (with respect to the principal stretch) of a strain energy function. Clearly, in this case 
the strain energy function must be dissipative and depend on the history of the strain, and so 
these authors express it as a fading-memory integral with integrand given as a hyperelastic strain 
energy function. This approach may be somewhat hard to justify and does not allow the user to 
consider the problem in terms of an auxiliary instantaneous measure of strain, i.e. the effective 
elastic stress term TT^ given above. 

Perhaps, a more subtle point to those mentioned above is the choice of TT^ in (1.1). In this 
paper, the reasonable assertion is made that TT^ must be zero whenever the deformation is zero. 
However, this point is often not recognized, especially for incompressible materials, as can be 
seen, for example, by inspection of the integrands (setting the stretch equal to unity) in the articles 
by [15] (see eqns (8) and (12) therein), [16] (eqns (9) and (13)) and [17] (eqn (9)). Note that this 
requirement is satisfied in the one-dimensional models in [8,9] but not in [10]. If TT^ does not 
vanish for zero deformation then, in general, it can be expected that the actual stress field TT will 
be non-zero prior to the application of the imposed load or stretch, i.e. the solution will be non- 
causal! However, for incompressible materials, an arbitrary pressure term (Lagrange multiplier) 
can be adjusted so that the stress field is indeed causal, but then it must take a specific form at 
later times. That is, whenever the material is deformed, and then returns back to zero, the stress TT 
will instantaneously return to zero too. This is clearly an unrealistic consequence of the particular 
choice of model, for, as Fung states [2, page 228]: 'the tensile stress at any time t is equal to the 
instantaneous [elastic] stress response . . . decreased by an amount depending on the past history'. 
This issue is examined further in §4fl. 

There are other variants of QLV employed in the literature that offer slightly modified forms of 
the governing equation to that suggested by Fung. The reader is referred, for example, to articles 
[18-20]. In this paper, the method employed follows closely that proposed by Fung but does not 
exhibit any of the limitations just discussed. 

The paper is organized as follows. In the following section, the usual (Boltzmann) linear 
viscoelastic model is reviewed by the way of introduction to a new interpretation of Fung's 
QLV theory, presented in §3. In §4, results from the new model are contrasted with those in 
the literature for the case of uniaxial elongation of a bar: in §4^ for an imposed stretch of an 
incompressible bar, in §4^ for an imposed load, and in §4c for stretch of a compressible bar. 
A numerical solution to a Volterra integral equation is required for the results in §4^, and the 
discretization (time-stepping) procedure used is described in appendix A. Finally, concluding 
remarks are offered in §5. 

2. Boltzmann's linear viscoelastic law 

It is helpful to commence analysis by recapping the theory of linear viscoelasticity. Under 
the assumption of isotropy infinitesimal elastic deformations can be described by the 
constitutive law 

u(t) = l^{e{t) - I tr(e(f))I) + k tr(e(0)I, (2.1) 

where cr and e are the second-order stress and strain tensors, respectively, and ji and k represent 
the infinitesimal shear modulus and modulus of compression (or bulk modulus), respectively. 



To incorporate viscoelastic behaviour, the most natural extension of (2.1) is to assume that the 
stress remembers the past history of the rate of strain with some fading memory, and then apply the 
superposition principle (or Boltzmann's principle); hence. 



(7(f) = 2 



M(f-s)— e(s) 
, 9s ' 



-(tre(s))ljds + 



K{t- 



■s)-(tre(s))Ids, 

9s 



(2.2) 



where /x(f), K{t) are now time-dependent relaxation functions, and the lower limit of the integral 
must be taken from — oo in order to correctly consider the initial deformation. Integrating (2.2) by 
parts, and assuming that the deformation commences at f = 0, yields 



(7(f) = 2M0)(e(0--(tre(0)1 1+2 



j + 2 IJi'it - s) {e{s) - ^(tr e(s))I j ds 



+ k{Q) tr e(f)I + 



;c'(f-s)tre(s)Ids, 



(2.3) 



where the ' denotes differentiation with respect to the argument of the function. Note that this 
expression incorporates any jump discontinuity when the motion starts. Moreover, the first term 
in (2.2) (or equivalently the first two terms in (2.3)) is trace free, or deviatoric, and accounts for 
shear deformations and loss, while the second term accounts for the hydrostatic part representing 
compressive deformations and loss. 

Now, in the elastic case (2.1), incompressibility may be considered as the dual limit of /c//x ^ oo 
and tr e ^ 0, which results in the stress-strain relationship 



u{t) = l^ie{t)-p{t)l, 



(2.4) 



i.e. the second term in (2.1) yields a finite non-zero limit, where p{t) may be considered 
as a Lagrange multiplier of incompressibility. For the viscoelastic counterpart, the limit of 
incompressibility can be considered in a similar fashion, noting first that /x(oo) and a:(oo) are the 
long time shear and bulk moduli (the equivalent of /x, k in the elastic case), respectively. Thus, 
taking a:(oo) -> oo (which implies /c(f) ^ oo for all t owing to the fading memory assumption) and 
tr e ^ 0 in (2.2) (or in (2.3)), it is found that 



k{0) tre(f) + 



K'{t - s) tr e(s) ds -p{t), 



(2.5) 



where the limit is again assumed to have a non-zero finite value, —p{t). Thus, in the limit of 
incompressibility, equations (2.2) and (2.3) become 



u{t) = -p{t)l^l 



f^{i-s) — e{s) ds 

) ds 



and 



cj{t) = -p{t)l + 2ix{0)e{t) + 2 



IJi'{t — s)e(s) ds. 



(2.6) 



(2.7) 



respectively. 



3. Quasi-linear viscoelasticity 

When the strain is not infinitesimal, linear theory becomes inappropriate to describe 
deformations; hence, a nonlinear constitutive law has to be considered. As discussed in the 
Introduction, Fung's hypothesis (1.1) is examined here as a means of describing the motion of 
viscous nonlinearly elastic materials. In his renowned work [2], Fung introduced this quasi-linear 
constitutive model in order to capture the nonlinear stress-strain relationship of living tissues; 
however, it also has applicability to elastomeric materials. 



Before deriving the quasi-linear theory, it is useful to introduce some standard definitions and 
notations. The deformation gradient tensor F is defined by 



F(S): 



I, se(-oo,0), 

ax(s) , , (3.1) 
se[0,f]. 



ax 

with x(s) denoting the position of a generic particle P at time s e [0, f], and X its position at the 
initial reference time. Note that the start time of the deformation, and any imposed tractions, will 
be taken as t = 0. The quantity / = det F, expressing the local volume change, is a constant J =1 
when the deformation is isochoric. Furthermore, from the deformation gradient tensor F the left 
Cauchy-Green tensor B = FF^ is obtained, together with its principal invariants 

Ji=trB, J2 = ^[(trB)^ -trB^] = (detB)tr(B-^) and l3 = detB=f, (3.2) 

which, alternatively, can be expressed in terms of the principal stretches through 

I^=xl-\-xl-\-kl, l2 = xlkl-\-klkl-\-klkl and Js^A^A^A^. (3.3) 

Now, Fung [2] makes the assumption that the QLV stress depends linearly on the (superposed) 
time history of a related nonlinear elastic response (a nonlinear instantaneous measure of strain). 
This allows, for example, for incorporation of a finite hyperelastic theory in the analysis. In index 
notation, Fung's theory (1.1) can be written as^ 



Gijkiit -t)^^ dr. (3.4) 

Fung refers to G/y^/ as a reduced relaxation function tensor. The 'crucial' simplification of Fung's 
theory is that this term is independent of the strain. Moreover, if the material is isotropic then G, a 
tensor of rank four, has just two independent components, being therefore consistent with linear 
theory.^ 

Following the analysis of the previous section, for isotropic materials it is convenient to split 
the equivalent (instantaneous) Cauchy stress into two parts, one which accounts for microscopic 
isochoric deformations of the material and the other that measures purely compressive deformations. 
These two components can be expected to have different instantaneous elastic behaviours as well 
as distinct relaxation rates, associated for example with the unwinding /unravelling of polymeric 
fibres as opposed to their stretching. It is assumed that this decomposition can be achieved by 
taking the deviatoric and hydrostatic components of the equivalent elastic Cauchy stress: 

T^ = T^ + Tfj, (3.5) 

which can be expressed as 

Tfj = ^ tr(T^)I and = - ^ tr(T^)I. (3.6) 

The split of equation (3.5) into shear and dilatational components is the nonlinear (hyperelastic) 
analogue of linear stress-strain law (2.1). However, it cannot be generalized to viscoelasticity, as 
in (2.2), because it would not preserve objectivity. Instead, the second Piola-Kirchhoff stress tensor 
associated with (3.5) has to be introduced, which is defined by 

= /F-^T^F-^ (3.7) 

with 

n^ = n[, + n^, (3.8) 

in which 

n^^/F-^T^F-T and ^/F-^T^F-^. (3.9) 

^Note, in order to clarify exposition, the space variable in (1.1) is omitted; the same will be done henceforth, as long as no 
confusion occurs. 

^This is not always taken in account, for example in [1,21] the authors introduce an isotropic model which contains only one 
scalar relaxation function G. 



It must be emphasized that the subscripts D and H do not refer to the deviatoric and hydrostatic 
parts of the second Piola-Kirchhoff stress, but correspond to the second Piola-Kirchhoff stress of 
the deviatoric and hydrostatic Cauchy stress components, respectively Assuming a superposition 
principle as for the linear case, it is now possible to introduce an objective viscoelastic law, relating 
the second Piola-Kirchhoff stress to the past history of the nonlinear rate of strain measure. This 
is taken as 



U{t)-. 



p(f-s)^n^(s)ds+ 



n{t-s)-ui,{s)ds, 



(3.10) 



where now T>(t) and H{t) are two scalar (independent)-reduced relaxation functions (with P(0) = 
H(0) = 1 without loss of generality). The latter relaxation functions relate to the inherent viscous 
processes involved with instantaneous shear and compressional (volumetric) deformations, 
respectively. Clearly, if the material was anisotropic, then a more complex tensorial relaxation 
function would be required. Furthermore, pre-multip lying by /~^F and post-multiplying by 
yield the Cauchy viscoelastic stress 

T(0 =r'F(f) D(f - s)^n^(s)ds) fT(0 
+ /-1f(0 (J^^ Hit - s)^n^(s) ds) fT(0, 

and integrating by parts, following (2.7), gives 

T{t)=J-'m (ul,{t) + V\t - s)n^(s) ds) ¥^{t) 



(3.11) 



+r'F(o (ui,{t) + \l^n\t - s)nfj(s)ds) i\t). 



(3.12) 



As mentioned earlier, to allow for large (nonlinear) deformations, a hyperelastic theory can 
be employed assuming that the measure of the effective elastic stress TT^ is derived from an elastic 
potential. Let us then specialize equations (3.5)-(3.12), considering the existence of a strain energy 
function (SEP) W, say, which in the isotropic case is dependent on the principal invariants of the 
deformation Ji, J2, 13, or of the principal stretches Ai, ^2/ -^3^ 



W = W(Ji, J2, 13) = W(Ai, A2, A3). 
The general form of elastic Cauchy stress may be written (e.g. [22,23]) as 

r = ^0l + iSiB + ^-iB-\ 



(3.13) 



(3.14) 



where = Pjih/h/h) are the so-called material (or elastic) response functions. In terms of the 
strain energy function, they are given by 



and 

where 



Mh,h,h) = jU2W2 + hW3l 
p.i{h,h,h) = -2m, 



(3.15) 



9W 



fc= 1,2,3. 



It is straightforward to calculate the trace of T*^, 



trr=3/Jo + Iift + 



J2/8-1 
/3 



(3.16) 



(3.17) 



and so the deviatoric and hydrostatic elastic Cauchy stress components in (3.6) become 



and 



= J ^lihm - hm)i + wiB - J3W2B-1 

rfj = y Q12W2 + ^II Wi + JsWs) I. 



,3 3 

From this, the second Piola-Kirchhoff counterparts, (3.9), are given by 



(3.18) 



and 



^(J2W2 - II Wi)C-^ + Wil - J3W2C-^j 

nfj = 2 Qj2 W2 + ^ Ji Wi + Js c-\ 



(3.19) 



(3.20) 



where it is recalled that C = F^F is the right Cauchy-Green deformation tensor. The viscoelastic 
stress is obtained from (3.19) to (3.20) via (3.10) or for the Cauchy stress from (3.12). 

Note that (as required) objectivity is now preserved for the viscoelastic model (3.11), but the 
first part is not deviatoric, and so the second term is not purely hydrostatic. In fact, in general 
both the compressive and shear components of the stress history contribute to the deviatoric and 
hydrostatic parts of T. However, the main point to note is that when there is no deformation, i.e. 
the principal stretches are unity, Ay = 1,; = 1, 2, 3 and B = I, then Ii =12 = 3,1^ = 1, which reveals 
that the effective deviatoric elastic stress vanishes. Similarly, T|j vanishes as the strain energy 
function has always to satisfy the additional relation (e.g. [22]) 



Wi(3, 3, 1) + 2W2(3, 3, 1) + W3(3, 3, 1) = 0. 



(3.21) 



Thus, both terms TT^ and TT^ vanish as 1, and hence, as there is no applied stress until the 
initial time t = 0, justifies taking the integration range for the pair of integrals in (3.12) as 0 to t. 

If it proves convenient to express the strain energy function in terms of the principal stretches 
W = W(Ai,A2,A3), then for diagonal F the quasi-linear viscoelastic stress relation (3.12) may be 
rewritten as 

Tiit) = r'xjit) (w^iit) + V'{t - s)n^,.(s) ds) 

0 / 



+}-'xHt)(uut) + 



(3.22) 



in which 



Di ■ 



A,- 



3Xf 



and 77' 



Hz " 



1 

3Xf 



7=1 



(3.23) 



where now W/ refers to the derivative dW/dXi. 

The final consideration of this section is the constraint of incompressibility for all possible 
deformations, i.e. /=1. In this limit, the speed of propagation of compressive disturbances 
tends to infinity, and hence the relaxation time for viscous dilatational motions can be 
assumed to tend to zero. Therefore, the second term in (3.11) (or (3.12)), in an analogous fashion 
to that for linear elasticity (2.7), reduces to 



j-'m (uf,{t) + n\t - s)nfj(s) ds) ¥^{t) ^ -p{t)i 



(3.24) 



where p{t) can be considered as a Lagrange multiplier. Note that this expression in itself is not 
hydrostatic. It contains both a hydrostatic term, and a component that combines with the first 



term in (3.12) to make it deviatoric; hence. 



T{t) = m (ul^it) + ^ V'{t - s)n^(s) ds^ I^{t) - pi (3.25) 
where now from the first equation in (3.19), 



n^(f) = 2 (^|w2-|wi^c-i + Wii-W2C- 



(3.26) 



As described above, the present theory has been developed for isotropic materials. The study of 
anisotropic nonlinear viscoelastic materials rapidly becomes rather demanding, although strain 
measures for anisotropic QLV were discussed in [24]. It is expected that, despite the complexity, 
the present approach is extendable to certain classes of anisotropy. 



4. Uniaxial elongation 



One of the simplest, useful and hence most popular experiments to measure the properties of a 
material is to subject a specimen to a simple elongation test. There is a huge literature of results on 
such a deformation, and it is therefore useful to examine this homogeneous solution, comparing 
the present result with extant published work. To appreciate the general QLV model developed 
herein, and the approach needed to obtain a solution, three specific cases will be examined. 
First, in §4fl a uniaxial stretch is imposed on the specimen, with the stress determined as a 
relaxing function of the history of the stretch. For comparison, both Yeoh and Mooney-Rivlin 
incompressible strain energy functions are considered. In §4^, a tensile load is imposed, and the 
stretch history (creep) determined from this. In this case, the integral equation must be solved 
numerically, which is the focus of discussion in appendix A. Finally, in §4c, the procedure in 
§4a is repeated for a compressible material, with the effect of compressibility on the viscoelastic 
behaviour highlighted. 

(a) Simple extension 

In the case of simple extension, assumed homogeneous (in space) and incompressible, the 
principal stretches may be specified as 

xi{t) = Xi{t)Xi, X2{t) = X2{t)X2 and X3(0 = ^2(0^3, (4.1) 

where (Xi,X2,X3) and (xi,X2/^3) are the Cartesian coordinates in the undeformed and deformed 
state, respectively, and the stresses are 

Tn{t) = m, T22{t) = T33{t) = 0 and Tij = 0{i^jl (4.2) 

having assumed that the lateral surfaces are free of stress. The deformation gradient Fij = dxi/dXj 
is of diagonal form 

F(f) = diag(Xi(f),X2(f),X2(0) (4.3) 



and the constraint of incompressibility, / = 1, gives a relationship between the stretches Xi and X2, 

— 1/2 

in particular X2 = '^i • Setting Xi=X, the deformation gradient tensor F and the left Cauchy- 
Green tensor B = FF^ become 

F(f) = dmg{x{t),X-^'^{t),X-^'^{t)^ and B(f) = dmg{X^{t),X-^{t),X-^{t)). (4.4) 

The principal invariants are therefore 

9 2 1 

II = A.^ + -, J2 = 2A + ^ and h = h (4.5) 

A A, 

and hence taking the diagonal terms of equation (3.25) yields the principal Cauchy stresses 

Tn(0 = m = X\t) (n^n(i) + V\t - s)n^,,{s) ds) - p (4.6) 



and 



Tiiit) = T33(f) = 0 = x-\t) ( n$,,Jt) + 



P'(f-s)/7^22(s)ds -p. 



(4.7) 



respectively. Thus, by subtraction, the Lagrange multiplier p can be eliminated, and so it is 
possible to rewrite the stress-strain relationship (4.6) and (4.7) as 

1 



T{t) = x\t)n^,,{t)---n^^^{t) 



where from (3.26) 



and 



Wi + 
Wi + 



W2 



1 

X3 



(1-/3) 



(4.8) 



(4.9) 



(4.10) 



Note that equations (4.9) and (4.10) depend on the specific choice of strain energy function W, 
and so it is useful to examine two specific examples. Let us start by assuming the instantaneous 
response is modelled by a two-term Yeoh strain energy function [25] 



W=|(2(/i-3)+a(/i-3)2), 



(4.11) 



which yields the stress-strain relationship 

T = -j^I + /x(l - 3a + aJi)B, (4.12) 
where a is a positive constant and /x is the shear modulus from infinitesimal theory. Thus, 



2Wi = fi{l - 3a + ah), W2 = 0. 
It is straightforward to obtain from (4.8) the relation 
T(f) 



(4.13) 



--m(m 



' A2(f)) 



D'(f - s)fc(s) 



where 



k{t) = 2a + (1 - 3a)X(t) + aX^{t). 



(4.15) 



The second example is the instantaneous response modelled by a Mooney-Rivlin strain energy 
function 



IX 1 



2 V2 



W=^U + y (/i-3)+^ 



2 \2 



y](h- 3), 



(4.16) 



which yields the stress-strain relationship 

T=-pI + ^[(l +y)B-(l-K)B-l], (4.17) 

where y is a constant in the range — 1/2 < y < 1/2. Note that (4.17) reduces to the neo-Hookean 
model when y = 1/2. The partial derivatives of W are 



2Wi=/i(l + y) and 2W2 = ii[^-y), 
and so (4.8) becomes, after simplification. 



(4.18) 




Figure 1. The imposed stretch history is shown in graph (a). The resultant dimensionless stress T/fi is plotted in graph (b), 
found for the Yeoh model (4.14) (dotted) and for the Mooney-Rivlin material (4.19) (dashed) where the solid curve is the 
neo-Hookean limits = 0 (or y =1/2). (c) T/fi is plotted from the predictions of (4.24) (dotted), (4.23) (dashed) and (4.14) 
(solid), respectively, (d) The dimensionless stress, T/fi, is plotted against stretch, X, from the predictions of (4.24) (dotted), 
(4.23) (dashed) and (4.14) (solid). 



in which 

m = [l-2y + X{t){l+2y)]. (4.20) 

These two viscoelastic models can be compared by considering an 'experiment' where the 
stretch is imposed and the stress is measured. To specify matters, the stress relaxation function is 
chosen to be the classical one-term Prony series 

V(t)=^^(l-^)e-'/\ (4.21) 

where /x, /Xoo are the infinitesimal shear modulus and the long-time infinitesimal shear modulus, 
respectively, and r is the relaxation time. These are set as 

— =0.5 and r = 1.0s. (4.22) 

A dynamic imposed stretch history (shown in figure la) is applied. The time variation is assumed 
slow so that inertial terms in the balance equations can be neglected. Figure lb shows the resultant 
stress predictions T//x for the Yeoh hyperelastic model (4.14) (dotted curves with a = 1,1) and for 
the Mooney-Rivlin hyperelastic model (4.19) (dashed curves with y = 1/6, —1/3). It is interesting 
to observe that the two hyperelastic models depart from the solid curve, obtained when the 
strain energy function is of the neo-Hookean type, i.e. a = 0, y = 1/2, in different ways. The Yeoh 
material is found to harden as the parameter a increases, whereas the effect of decreasing y from 
1/2 leads to a softening of the Mooney-Rivlin material's behaviour. 

As mentioned in the Introduction, it is useful to contrast the results presented here with those 
discussed, for example, by Ciambella et al. [15]. In that article, the authors employed a form of 



QLV, which, with the incompressible Yeoh SEF, yields the stretch-stress law (see eqn (16) in [15]): 
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— = [k{t) - X-^{t)][2a + (1 - 3a)k{t) + ak^{t)] 

rt 

+ [X^{t) - k-^{t)] V'{t - s)X-^{s){la + (1 - ?>a)X{s) + aX^{s)] ds. (4.23) 
Jo 

Similar equations can be derived from the constitutive models in [16] or [17]. In Ciambella et al.'s 
article, their equation (4.23) was compared with that offered by the formulation used in the 
ABAQUS finite-element analysis (FEA) package (Version 6.7, see [26]), namely 

^ = [X(t) - k-^(t)]l2a + (1 - 3a)k(t) + ak^(t)] : J 

t 

V'{t - s)[A(s) - X-^{s)]{la + (1 - 3oi)X{s) + aX^{s)] ds. (4.24) 

0 

The three alternative viscoelastic expressions (4.24), (4.23) and the present result (4.14) may be 
conveniently compared by again imposing a stretch and determining the resultant dimensionless 
stress T//X. Employing the same stretch history as above (figure \a), the relaxation function in 
(4.21) with parameters as in (4.22) and a = 1.0, yields the curves in figure Ic predicted from 
(4.24) (ABAQUS: dotted), (4.23) (Ciambella et al\ dashed) and (4.14) (De Pascalis et al: solid), 
respectively. An alternative way of representing this deformation is via a stretch-stress diagram 
(figure Id). It is clear that the present QLV approach, illustrated by the solid curve in both 
figures, gives a result that lies remarkably close to that found from the ABAQUS model, whereas 
Ciambella et al.'s solution is quite distinct. The latter model was developed by the authors in order 
to address the deficiency they highlighted regarding v. 6.7 of ABAQUS. However, their result also 
exhibits a substantial limitation: it predicts instantaneous zero stress whenever X recovers back to 
unity after some (arbitrary) deformation, thus showing no 'memory' of the history of the stress in 
this situation. This is at variance with physically observed behaviour, and in particular that found 
for linear viscoelasticity. 



<=> 
<=> 

on 
oo 



(b) Simple tensile load 

Rather than an imposed stretch, perhaps a more important experiment for measuring viscoelastic 
properties of materials is application of a simple tensile load: a slowly varying uniaxial stress 
is imposed on the body for which the resultant stretch is measured over time. For linear 
viscoelasticity, equation (2.2) can be inverted to yield a straightforward creep relation to determine 
the strain. However, for nonlinear viscoelasticity, the solution procedure is somewhat more 
difficult, as inversion is not possible. A typical QLV relation is given in (4.19), which may be 
considered in the general form 



.0 



m =gm) + Y^fjm) ^ T^\t - s)hj{x(s)) ds. (4.25) 



Although the integrand is of separable type, the presence of the various ^- (A (f)) terms means that 
an inversion operator cannot be introduced. Instead, a numerical scheme must be employed that 
can evaluate the stretch as a function of time, subjected to a prescribed stress. A suitable numerical 
discretization procedure is offered in appendix A, which has error 0{8t^), where 8t is the step size 
and so gives rapid convergence. 

An example is illustrated in figure 2, with the imposed stress history shown on the left graph. 
The one-term Prony series relaxation function (4.21), with constant values as given in equation 
(4.22), is again employed. The resultant stretch is given in figure 2h for the Yeoh strain energy 
function (4.14) (dotted curves with a = 1,2), and for the Mooney-Rivlin strain energy function 
(dashed curves with a = 1,2). As before, the solid curve is the neo-Hookean prediction {a = 0, 



0 5 10 15 20 25 30 0 5 10 15 20 25 30 

t t 



Figure 2. Plot of the dimensionless stress history T/ii [a) and the resultant stretch A. [b), from the Yeoh model predictions 
(4.14) (dotted) and that from Mooney-Rivlin predictions (4.19) (dashed). The solid curve is the neo-Hookean limit a = 0 
(orK=1/2). 



y = 1/2) and it can be seen that increasing a in the Yeoh model leads to material hardening, while 
decreasing / leads to softening of the Mooney-Rivlin material. 



(c) Simple extension of a compressible bar 

The purpose of this section is to repeat that analysis in §4fl for a compressible material. As before, a 
bar, with zero tractions on its lateral faces, is assumed to undergo simple homogeneous extension, 
caused by an imposed stress or stretch in the Xi direction, say. Then, the principal stretches may 
be given by 

xi{t) = Xi{t)Xi, X2{t) = k2{t)X2 and xsW^ ^2(0^3, (4.26) 

where symmetry between the X2 and X3 directions is clear, and so the deformation gradient 
tensor is 

F(f) = diag(Ai (0,^2(0,^2(0) (4.27) 



and the left Cauchy-Green tensor B becomes 

B{t)^diag{xl{tlxl{t)Alm (4.28) 

The principal invariants are therefore 

h{t) = xl{t) + 2xl{t\ I2{t) = 2xl{t)xl{t) + Xj{t) and h = xl{t)xj{t). (4.29) 

The above can be substituted into the governing viscoelastic equation (3.12), with effective 
elastic stresses (3.19) and (3.20), and simplified. The Xi component of the stress yields, after 
simplification. 



T(0 = 



Ai(0 

xl(t) 



D'(f - s) 1 



A|(s) 
Xlis) 



(Wi(s) + W2(s)X^(s))ds 



+ 



+ 2 



W'(f-s) Wi(s) + 2 



mis) 



Wi(s) 
Xl(s) 



+ 2W2(s)U^(s) 



xjis) 



+ 3W3(s) A*(s) ds + 2(Wi(0 + 2W2(0^2(0 + ^3(0^2(0) 



(4.30) 



and the normal stresses in the X2 and X3 equations are both 



0 = 



V'(t - s) 



0 \ ^2(S) 



+ 2(Wi(s) + 2W2(s)Af(s)) 



+ (2W2(s) + 3W3(s)a2(s))a2(s) ) ds 



+ 2(Wi(f) + W2(f)^?(f) + (W2(0 + W3(0A?(f))^2(0) 



(4.31) 



where as before Wy is the derivative of W with respect to the indicated invariant (3.16). Clearly, 
the particular choice of the strain energy function affects the results and so, as before, a couple 
of examples are presented. Assuming the Horgan-Murphy strain energy function [27], used for 
modelling a material with a small compressibility, the SEF is 



W-- 



32) 



where y is an arbitrary constant, /x is the infinitesimal shear modulus and k is the infinitesimal 
bulk modulus. Then 



and 



^r'^'-Mf^-K)/3-^^v|(i-/3-n 



which are substituted in (4.30) and (4.31). Setting y = 1/2, for ease of presentation, the stresses 
become 



T(f) = 



Kit) IX + 



2/x 



I ^?(s) 



•a2(s) 



ds 



n'{t - s) ( M ( 2^ + 1 - 3^1^ 

(s) 



+ 3k I x|(s) - :jl5fl ) ) ds + ,fx|(f) ) - xl'\t)x'^2^^(t)(ix + Kxy\t)xl'^(t)) 



his) 



(4.33) 



and 



0 = 



Xi(t)Xl(t) 



^2(0 ( + f 



D'(t-s)ll--A^ Ids 
0 \ A|(s) ^ 



7i'(f - s) /X 2 + 



A^(s) X^/3(s) 



+ 3/f Ai(s)(Ai(s)A5(s) - 1) ds 



+ Kxl(t)xl(t) - ^f)^2^^0(M + Kxy^(t)xl^ht)) 

respectively. By contrast, for the Gent model (e.g. [28]), 



(4.34) 



where Jm is a constant limiting value for Ii — 3, taking into account the finite extensibility of 
polymeric chains within the material. Hence, 



which can then substituted into (4.30) and (4.31) to yield the equivalent results for the stresses. 
These are omitted for brevity. 



This paper has focused on reappraising Fung's method for QLV. It has been shown that some 
of the negative features commonly associated with the approach are merely a consequence of 
the way it has been applied elsewhere. The approach outlined herein was shown in §4 to yield 
'sensible' results and offers a straightforward approach in solving a wide range of models. 
The present method exhibits similarities with Simo's approach to nonlinear viscoelasticity [7], 
although the latter assumed a scalar relaxation function acting on the stress. Therefore, it is 
noted that Simo's method would not reduce in the linear limit to the usual relation (2.3) if the 
latter shear and bulk moduli have different temporal behaviour. (Note that ABAQUS v. 6.12 
employs the same relaxation time constants in both the shear and bulk parts of the field.) Herein, 
the relaxation function is a tensor, which for isotropic materials reduces to two distinct scalar 
relaxation functions, one acting on the compressive part and the other on the shear component of 
the stress. This therefore is consistent with that found for linear theory. 

It was shown that for an imposed deformation (e.g. stretch) it is simple matter to solve the 
QLV equation directly. For imposed stress, the stretch has to be deduced by solving (inverting) 
the integral equation, and this is achieved here using a discretized scheme accurate to an order 
of the cube of the time-step size. Muliana et al. [29] have recently offered a QLV model where 
the strain is expressed as a function of the stress, which may be viewed as a dual model to (1.1). 
However, it is not clear as yet whether their approach is as effective at modelling viscoelasticity 
as Fung's scheme. 

The authors are currently applying the present approach to a range of problems in 
viscoelasticity. The numerical procedure described in appendix A can be employed in a 
straightforward manner for any deformation that is incompressible and homogeneous, for 
example simple shear. It is presently being extended so that it can used for compressible materials 
undergoing simple extension or shear. The QLV method is also adaptable to inhomogeneous 
problems, for example those studied in [30], large amplitude radial deformations in two 
and three dimensions, and simple torsion. Such models lead to partial differential nonlinear 
Volterra integral equations and extending the present numerical scheme to these should offer 
improvement over current methods. It is also anticipated that the present approach will be 
more suited to satisfying the equations of equilibrium than, say, the more heuristic QLV relation 
employed in ABAQUS (v. 6.12). Finally, the authors are using the model for several biomechanics 
problems, to derive a perturbation theory for the viscoelastic evolution of small deformations on 
the top of large ones and to derive effective properties of ligaments and other tissues. 

Acknowledgement. The authors are grateful to the Engineering and Physical Research Council for the award 
(grant no. EP/H050779/1) of a postdoctoral research assistantship for R.DeP. 



Appendix A, A numerical high-order solution procedure for Volterra integral 
equations 



The purpose of this appendix is to offer a simple and rapid numerical scheme for solving Volterra 
integral equations of the type obtained by the QLV approach (see e.g. (4.14)). It is assumed that 




(4.35) 



5. Concluding remarks 



they have the separable form 



N t 

m =8{X{t)) + TfjiHt)) V'{t - s)hj{X{s)) ds, 



(Al) 



where T{t),V{t),g{X),fj{X) and hj{X) are all known functions of their respective arguments. 
Clearly, all the equations offered in §Aa lie within this class but for more general problems, with 
inhomogeneous deformations, the analysis described below must be amended. 

Let us now discretize the system and evaluate the equation at each time tn = nSt, with 
n = 0, 1, 2, 3, . . . , nmax/ where 8t is a small time step. On writing the stretch at each time step as 



kn=Htn) = Hn8t\ 

a Taylor series expansion can be employed to yield 

SiK) =g{Xn-l)+8'{K-l){X„ - X„-l) + ^ ^^~^' {Xn - X„-lf 
■ ^"'^^"-^\xn - Xn-lf + 0{{Xn - Xn-lf), 



(A 2) 



3! 

where ' denotes differentiation with respect to the argument, and 

X" X'" 

X„ = X{n&t) = A„_i + X'„_^&t + ^5*2 + ^5f3 + 0{&t\ 



(A3) 



(A 4) 



Note that all functions are assumed sufficiently smooth between adjacent time steps so that the 
relations above and below hold. From (A 4), it is clear that 



{Xn - = k^^_,8t + -^8f + -^8t^ + 0{8t\ 

and {Xn - K-lf = K-lf^i^ + 0{8t^)- 

The integral in (A 1) can also be expanded as 

t rnSt 

V'{t - s)hAX{s)) ds = V'{n8t - s)h{X{s)) ds 
0 Jo 



(A 5) 



n-l 
m=0 



mSt 



V'{n8t - s)hj{X{s)) ds 



= ^ P'((n - m)8t - u)hj{X{m8t + u)) du, 



m=0 



(A 6) 



where in the latter the substitution s = m8t + u is used. The integrand functions can be further 
expanded as 



V\{n - m)8t -u) = V\{n - m)8t) - V'\{n - m)8t)u 
, V"'{{n-m)8t), , 



2! 



(A 7) 



and 



<^hAXm) 1 Sh:{Xm) , , 

hj{X{m&t + u)) = hjiXm) + -^^u + - + 0{u^ 



-.hj{X,n) + X'Jl'-{X,n)u + -{X';„h'j(X,n) + {X'J^ h'j {X,n)V + 0{u^), 



(A 8) 



and therefore, after some simple algebra and integration, it is found that 

rSt 

V'{{n - m)8t - u)hAk{m8t + u))du = V'{{n - m)8t)hAkm)8t 

.0 

8t^ 

+ (V\{n - m)8i)y!j{.{\^) - V"{{n - m)8t)hj{km))^ 
+ {V'\{n - m)8t)hj{Xm) - 2V'{(n - m)8t)X'Ji'.{Xm) 

8t^ 

+ V\{n-m)8t)(Ch^.{krn) + {k^J^hJ^ (A 9) 

; b 

Using now the expansion given in (A 3), with the help of relation (A 5), it is possible to rewrite • 
the expansion for both g and fj as 

g{Xn)=g{Xn-l) + X'^_^g\kn-l)8t + {k';_^g\kn-l) + g'' {Xn-l)) 



8f- ' ^ 



8t^ 



and 



-1^ 
o 

on 
oo 



+ {K'_Jl{Xn-i) + 3X;_iA;'_i^"(A„_i) + (X'„_,ffl"{Xn-l))-^ + 0(Sf4) 



Finally, gathering all terms, whole equation (A 1) can be written as the following expansion: 

8fi 

T{n8t)=g(kn-l) + k'^_^g\kn-l)8t + {k';_^g\kn-l) + i^-lf g' (kn-l)) — 

8t^ 

+ {k^;^_,g\kn-i) + 3a;_ia;^_i/(a^_i) + (a;_,)V^(a^_i))— 

8t^ 8fi A 

+ An-i8t + — + Cn-i — + 0{8t^\ (A 10) 
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where C„_i are the following expressions: 

N n-1 

= E J2^M^n-i)V\{n - m)8t)hj{km)l (A 11) 

y=l m=0 
N n-1 

= E J2^K-^j(^n-i)V\{n - m)8t)hj{km) 

j=l m=0 

+ fj{kn-i){V\{n - m)8t)k'Ji'.{km) - V\{n - m)8t)hj{km))\ (A 12) 

N n-1 

^n-1 = E Y.^fj{>^n-i){V'"{{n - m)8t)hj{km) - lG"{{n - m)8t)k'Ji'.{km) 
j=l m=0 

+ V'{{n - m)St)iX';„h'j(X„) + {X'^fh'! {X,n))) 
+ 3X'„_Jl{Xn-i){V'{(n - m)St)X'^h'j{X,n) - V"({n - m)St)hj(X„,)) 
+ 3V'i(n - m)St)hj(Xm)(Xl_j!(X„_i) + iX'„_if ff (X„_i))]. (A 13) 



and 



= 177^ ^ • (A 16) 



This equation must hold for each n and at each level of approximation. Hence, keeping terms up 
to 0(<5f) in (A 10), and rearranging, gives 

and then equating terms at higher order (0(<5f^) and 0{8t^)) to zero yields, respectively, 
and 

(K-iK-is"(^n-i) + (a;_i )V"(a„-i) + c„_i) 

Now, the procedure for solving the discretized equations (A 10), for specified stress, is clear. 
Deformation commences at f = 0 (or w = 0) and so from (A 1): 

r(0)=g(A(0)) (A 17) 

gives 

X{0) = Xo=g-\W))- (A 18) 

As there is no stretch prior to an applied stress, an initial zero stress, T(0) = 0, would mean that 
A-o = 1. The derivatives of are derived from (A 14) to (A 16) (setting n = 1 in these equations) 
and hence Ai is found from (A 4). Then n is incremented, the values of X' , X" , X'" determined, and 
these are again used to determine the stretch at the next time step. The process is repeated at 
each following step. Note that size of the sums in An-i, Q-i grows with n, and so for large 
times, or using small time steps, evaluation of the system can become rather slow. However, this 
difficulty can be alleviated, with the present system, by obtaining a relationship between An and 
its previous increment and similarly for Bn,Cn. The precise forms of these relationships 

are determined from the given /J- (A) and the relaxation function V{t). In practice, the solution 
procedure works well, with an error of 0(8t^), i.e. for a time step of 8t = 0.01 then the error is 
0(10-8). 
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